/*
 * =============================================================
 * h2tc.c - 
 *
 *      Computes temperature (T) for specified enthalpy (H) 
 *               & fuel-air ratio (fa)
 *
 *        D K Frederick
 * =============================================================
 */

#include <stdio.h>
#include <math.h>

double h2tc(double H, double fa)
{

	/*----- define things we will need ------*/
	int ITAB[60] = {1,1,1,2,3,3,4,4,4,5,5,5,5,6,6,6,6,6,7,7,7,7,7,7,
		 8,8,8,8,8,8,9,9,9,9,9,9,10,10,10,10,10,10,10,10,
		 11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11};
	int TTTAB[11] = {300,400,500,700,1000,1400,1900,2500,3100,3700,4500};
	double AHAIR[11] = {2.074402000000e3,2.767580000000e3,3.461230000000e3,
		 4.854285000000e3,6.981156000000e3,9.933801000000e3,
         1.382117000000e4,1.870811000000e4,2.375618000000e4,
         2.891518000000e4,3.591874000000e4};
	double BHAIR[11] = {6.930375143000e0,6.933262304370e0,6.941415639521e0,
         6.999806554132e0,7.201291270059e0,7.564016167549e0,
         7.965574033453e0,8.298411651745e0,8.515829359567e0,
         8.673620909986e0,8.825346607306e0};
	double CHAIR[11] = {1.327409630363e-5,1.559751739222e-5,6.593583412638e-5,
         2.260187389251e-4,4.455969808316e-4,4.612152628951e-4,
         3.419004689130e-4,2.128288949057e-4,1.495339514650e-4,
         1.134519659006e-4,7.620515574828e-5};
	double DHAIR[11] = {7.744736961966e-9,1.677943891139e-7,2.668048413312e-7,
         2.439758243406e-7,1.301523505289e-8,-7.954319598805e-8,
         -7.170643000405e-8,-3.51638574671e-8,-2.004554753575e-8,
         -1.551950423014e-8,2.007781800899e-22};
	double APAIR[11] = {4.229854000000e1,4.429218000000e1,4.584067000000e1,
         4.818303000000e1,5.070949000000e1,5.318950000000e1,
         5.556056000000e1,5.779398000000e1,5.960312000000e1,
         6.112402000000e1,6.283719000000e1};
	double BPAIR[11] = {2.305525000000e-2,1.732687898835e-2,1.390113404660e-2,
         9.984237743689e-3,7.194810211645e-3,5.398110354070e-3,
         4.191444392129e-3,3.318644249749e-3,2.746778608874e-3,
         2.344441314755e-3,1.960623219312e-3};
	double AHSTOC[11] = {2.116286000000e3,2.831822000000e3,3.556413000000e3,
         5.033766000000e3,7.326000000000e3,1.054811000000e4,
         1.483580000000e4,2.028520000000e4,2.596610000000e4,
         3.180568000000e4,3.976139000000e4};
	double BHSTOC[11] = {7.113538900000e0,7.199470008152e0,7.292391067390e0,
         7.482468579354e0,7.811853002145e0,8.297006217518e0,
         8.832077018487e0,9.300821657633e0,9.616136350979e0,
         9.837032938452e0,1.003677698592e1};
	double CHSTOC[11] = {3.953219184767e-4,4.639891630481e-4,4.652214293286e-4,
         4.851661304903e-4,6.127819454796e-4,6.001010929530e-4,
         4.700405089856e-4,3.112005562576e-4,2.143239326515e-4,
         1.538370464703e-4,9.584301286375e-5};
	double DHSTOC[11] = {2.288908152379e-7,4.107554268465e-9,3.324116860284e-8,
         1.417953499881e-7,-1.056737710548e-8,-8.670705597826e-8,
         -8.824441818226e-8,-5.382034644785e-8,-3.360382565619e-8,
         -2.416418066940e-8,6.080973759837e-15};
	/*-----------------------------------------*/
	
	double zmea, zmsp, tmlsr, zmwtr, tmls, tg, tgo, hgo, z, zz, hg, hh, hgsp, d, hgea, d1, temp;
	int ii, Tindex, Tindex2, it;

	if (fa == 0){
		zmea = 4.7642;
		zmsp = 0;
		tmlsr = 0.2098988288;
		zmwtr = 0.0345194683;
	}
	else{
		zmea = 4.7642-fa*69.69056873;
		zmsp = fa*74.411931335;
		tmls = 4.7642+fa*4.721362582;
		zmwtr = tmls/(138.014872100*(1+fa));
		tmlsr = 1/tmls;
	}

	tg = 3.55*H;
	tgo = 1;
	hgo = 0.282;
	z = H*2e-5;
	zz = tmlsr*zmwtr;
	
	ii=0;
	hh=99.99;

	while(fabs(hh)>1e-3 && ii<11){
		ii = ii+1;
		temp = 0.01 * (tg - fmod(tg, 100));
		if (temp > 1)
			Tindex2 = temp;
		else
			Tindex2 = 1;
		if (Tindex2 < 60)
			Tindex = Tindex2;
		else
			Tindex = 60;

		it = ITAB[Tindex - 1] - 1;
		d1 = tg - TTTAB[it];
		hgea = ((DHAIR[it]*d1 + CHAIR[it])*d1 + BHAIR[it])*d1 + AHAIR[it];

		if (fa > 0)
			hgsp = ((DHSTOC[it]*d1 + CHSTOC[it])*d1 + BHSTOC[it])*d1 + AHSTOC[it];
		else
			hgsp = 0;
		hg = (hgsp*zmsp+hgea*zmea)*zz;
		hh = H - hg;

		d = (tg - tgo)/(hg - hgo);
		tgo = tg;
		hgo = hg;
		tg = tg + d*hh;
	}

	return tg;
}

